# ------------------------------------------------------------------------------#
# Description: 
# Create:
#   Table 5: Peer effects by type of technology adopted
#   Table 6: Peer effects by type of technology adopted and peer technology
# Project:     Peer Effects in Voluntary Environmental Policies:
#              An Application to Urban Water Quality
# Author:      Daniel A. Brent | dab320@psu.edu
# Created:     2025-06-11
# ------------------------------------------------------------------------------#

# Clear environment ------------------------------------------------------------
rm(list = ls())
gc()

# Load required packages -------------------------------------------------------
library(pacman)
p_load(data.table, fixest, modelsummary, stringr, lfe, stargazer)

options(scipen = 999)

# Load adoption data -----------------------------------------------------------
load("data/adoptions/adoptions.RData")

# Merge peer adoption buckets --------------------------------------------------
load("data/bucket/peer_adopt_buckets_any.RData")
dat <- merge(adopt, peer.adopt.bucket.any, by = c("pin", "year"))
rm(adopt, peer.adopt.bucket.any)

load("data/bucket/peer_adopt_buckets_rg.RData")
dat <- merge(dat, peer.adopt.bucket.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.rg)

load("data/bucket/peer_adopt_buckets_any_rg.RData")
dat <- merge(dat, peer.adopt.bucket.any.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg)

load("data/bucket/peer_adopt_buckets_cs.RData")
dat <- merge(dat, peer.adopt.bucket.cs, by = c("pin", "year"))
rm(peer.adopt.bucket.cs)

load("data/bucket/peer_adopt_buckets_both.RData")
dat <- merge(dat, peer.adopt.bucket.both, by = c("pin", "year"))
rm(peer.adopt.bucket.both)

# Merge peer eligibility buckets -----------------------------------------------
load("data/bucket/peer_elig_buckets_cs.RData")
dat <- merge(dat, peer.buckets.cs, by = c("pin", "year"))
rm(peer.buckets.cs)

load("data/bucket/peer_elig_buckets_rg.RData")
dat <- merge(dat, peer.buckets.rg, by = c("pin", "year"))
rm(peer.buckets.rg)

# Merge total peers ------------------------------------------------------------
load("data/bucket/total_peers.RData")
dat <- merge(dat, total.peers, by = "pin")
rm(total.peers)

# Merge parcel data ------------------------------------------------------------
load("data/parcel_build/parcel_build.RData")
dat <- merge(dat, parcel.build[, .(pin, zip5, sub.area, sqft, lot, beds, baths,
                                   yearbuilt, val.land, val.improve,
                                   ren.pre90, ren.90.99, ren.00.09, ren.10,
                                   water.prob)], by = "pin")
rm(parcel.build)

# Merge census data ------------------------------------------------------------
load("data/census/all_parcels_census_2010.RData")
parcels.2010[, blk := as.numeric(paste0(block, blkgrp))]
dat <- merge(dat, parcels.2010[, .(pin, tract, blkgrp, block, blk)])
rm(parcels.2010)

# Scale peer adoption variables ------------------------------------------------
dat[, peer.adopt.any.100.scale    := peer.adopt.any.100 / 100]
dat[, peer.adopt.any.rg.100.scale := peer.adopt.any.rg.100 / 100]
dat[, peer.adopt.rg.100.scale     := peer.adopt.rg.100 / 100]
dat[, peer.adopt.cs.100.scale     := peer.adopt.cs.100 / 100]
dat[, peer.adopt.both.100.scale   := peer.adopt.both.100 / 100]

# Keep post-2010 only ----------------------------------------------------------
dat <- dat[year > 2010]

# Create average eligibility variable ------------------------------------------
dat[, elig.peers.avg.100 := elig.peers.cs.100 + elig.peers.rg.100]

# Table 5: Peer effects by type of technology adopted ------------------------------------------------------------------------------

base <- feols(adopt.any ~ peers.100 | year + blkgrp | 
                peer.adopt.any.100.scale ~ elig.peers.avg.100,
              cluster = 'blkgrp', data = dat[elig.cs == 1 & post.rainwise == 0])

m1 <- feols(adopt.any.rg ~ peers.100 | year + blkgrp | 
              peer.adopt.any.100.scale ~ elig.peers.avg.100,
            cluster = 'blkgrp', data = dat[elig.cs == 1 & post.rainwise == 0])

m2 <- feols(adopt.cs ~ peers.100 | year + blkgrp | 
              peer.adopt.any.100.scale ~ elig.peers.avg.100,
            cluster = 'blkgrp', data = dat[elig.cs == 1 & post.rainwise == 0])

m3 <- feols(adopt.rg ~ peers.100 | year + blkgrp | 
              peer.adopt.any.100.scale ~ elig.peers.avg.100,
            cluster = 'blkgrp', data = dat[elig.cs == 1 & post.rainwise == 0])

m4 <- feols(adopt.both ~ peers.100 | year + blkgrp | 
              peer.adopt.any.100.scale ~ elig.peers.avg.100,
            cluster = 'blkgrp', data = dat[elig.cs == 1 & post.rainwise == 0])

# Format IV diagnostics --------------------------------------------------------
load("data/iv/iv_diag_table3_4.RData")
eff.f.1 <- round(iv_diag_table3_4[["table3.iv"]]$F_stat[4], 2)
eff.f.2 <- round(iv_diag_table3_4[["table4.any.rg"]]$F_stat[4], 2)
eff.f.3 <- round(iv_diag_table3_4[["table4.cs"]]$F_stat[4], 2)
eff.f.4 <- round(iv_diag_table3_4[["table4.rg"]]$F_stat[4], 2)
eff.f.5 <- round(iv_diag_table3_4[["table4.both"]]$F_stat[4], 2)

ar.ci.1 <- paste0("[", round(iv_diag_table3_4[["table3.iv"]]$AR$ci[1], 2), ", ",
                  round(iv_diag_table3_4[["table3.iv"]]$AR$ci[2], 2), "]")
ar.ci.2 <- paste0("[", round(iv_diag_table3_4[["table4.any.rg"]]$AR$ci[1], 2), ", ",
                  round(iv_diag_table3_4[["table4.any.rg"]]$AR$ci[2], 2), "]")
ar.ci.3 <- paste0("[", round(iv_diag_table3_4[["table4.cs"]]$AR$ci[1], 2), ", ",
                  round(iv_diag_table3_4[["table4.cs"]]$AR$ci[2], 2), "]")
ar.ci.4 <- paste0("[", round(iv_diag_table3_4[["table4.rg"]]$AR$ci[1], 2), ", ",
                  round(iv_diag_table3_4[["table4.rg"]]$AR$ci[2], 2), "]")
ar.ci.5 <- paste0("[", round(iv_diag_table3_4[["table4.both"]]$AR$ci[1], 2), ", ",
                  round(iv_diag_table3_4[["table4.both"]]$AR$ci[2], 2), "]")

varnames <- c('year' ='Year','blkgrp' = 'Block Group',
              'peer.adopt.any.100.ols'= 'Peer Adoptions',
              'elig.peers.cs.100' = '\\# Eligible Peers (CS)',
              'elig.peers.rg.100' = '\\# Eligible Peers (RG)',
              'peer.adopt.any.100'= '$\\widehat{\\text{Peer Adoptions}}$',
              'peer.adopt.any.100.scale'= '$\\widehat{\\text{Peer Adoptions}}$')

# Preview Table --------------------------------------------------------------
etable(base, m1, m2, m3, m4,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE,
       dict = varnames,
       signif.code = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
       headers = list("_ " = c("Any", "Rain Garden", "Only Cistern", "Only Rain Garden", "Both")),
       fitstat = c("n"),
       extralines = list(
         "_F"     = c(eff.f.1, eff.f.2, eff.f.3, eff.f.4, eff.f.5),
         "_AR CI" = c(ar.ci.1, ar.ci.2, ar.ci.3, ar.ci.4, ar.ci.5)
       ))

# Save Table 5 --------------------------------------------------------
etable(base, m1, m2, m3, m4, tex = TRUE,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE,
       dict = varnames,
       signif.code = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
       headers = list("_ " = c("Any", "Rain Garden", "Only Cistern", "Only Rain Garden", "Both")),
       fitstat = c("n"),
       extralines = list(
         "_F"     = c(eff.f.1, eff.f.2, eff.f.3, eff.f.4, eff.f.5),
         "_AR CI" = c(ar.ci.1, ar.ci.2, ar.ci.3, ar.ci.4, ar.ci.5)
       ),
       file = "output/tables/table_5.tex", replace = TRUE)

rm(base, m1, m2, m3, m4, iv_diag_table3_4,
   eff.f.1, eff.f.2, eff.f.3, eff.f.4, eff.f.5,
   ar.ci.1, ar.ci.2, ar.ci.3, ar.ci.4, ar.ci.5)

# Table 6: Peer effects by type of peer technology ------------------------------------------------------------------------------

m1 <- felm(adopt.any ~ peers.100 | blkgrp + year |
             (peer.adopt.any.rg.100.scale | peer.adopt.cs.100.scale ~ elig.peers.cs.100 + elig.peers.rg.100) | blkgrp,
           data = dat[elig.cs == 1 & post.rainwise == 0])

m2 <- felm(adopt.any.rg ~ peers.100 | blkgrp + year |
             (peer.adopt.any.rg.100.scale | peer.adopt.cs.100.scale ~ elig.peers.cs.100 + elig.peers.rg.100) | blkgrp,
           data = dat[elig.cs == 1 & post.rainwise == 0])

m3 <- felm(adopt.cs ~ peers.100 | blkgrp + year |
             (peer.adopt.any.rg.100.scale | peer.adopt.cs.100.scale ~ elig.peers.cs.100 + elig.peers.rg.100) | blkgrp,
           data = dat[elig.cs == 1 & post.rainwise == 0])

m4 <- felm(adopt.rg ~ peers.100 | blkgrp + year |
             (peer.adopt.any.rg.100.scale | peer.adopt.cs.100.scale ~ elig.peers.cs.100 + elig.peers.rg.100) | blkgrp,
           data = dat[elig.cs == 1 & post.rainwise == 0])

m5 <- felm(adopt.both ~ peers.100 | blkgrp + year |
             (peer.adopt.any.rg.100.scale | peer.adopt.cs.100.scale ~ elig.peers.cs.100 + elig.peers.rg.100) | blkgrp,
           data = dat[elig.cs == 1 & post.rainwise == 0])

# Preview Table --------------------------------------------------------------
stargazer(m1, m2, m3, m4, m5, type = "text", keep = c("peer.adopt"))

# Save Table 6  --------------------------------------------------------
stargazer(m1, m2, m3, m4, m5, type = "latex",
          out = "output/tables/table_6.tex",
          float = FALSE, omit.table.layout = "n",
          dep.var.caption = "", dep.var.labels.include = FALSE, model.names = FALSE,
          keep = c("peer.adopt"), omit.stat = c("f", "ser", "adj.rsq", "rsq"), no.space = TRUE,
          column.labels = c("Any", "Rain Garden", "Only Cistern", "Only Rain Garden", "Both"),
          covariate.labels = c(
            "$\\widehat{\\text{\\ Peer Rain Garden Adoptions }}$",
            "$\\widehat{\\text{\\ Peer Cistern Adoptions}}$"
          ))
